{
 "metadata": {
  "name": "",
  "signature": "sha256:828ac6583c93644eaac78c8b02a0eb5403e615d12d56d41b201fc952b0c4780d"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Introduction to finite element method\n",
      "\n",
      "##Introduction\n",
      "\n",
      "The *Finite element method* (FEM) is a numerical procedure to obtain approximate solutions of boundary value problems governed by a differential equation. It basically produces approximate values of the desired parameters at certain points called nodes. \n",
      "\n",
      "The main characteristics that distinguish it from other numerical technics are:\n",
      "\n",
      "1. Utilizes an integral formulation to generate a system of algebraic equations.\n",
      "\n",
      "2. Uses continuous piecewise smooth functions to approximate the unknown quantities.\n",
      "\n",
      "The differences between FEM and FDM are: \n",
      "\n",
      "* FEM is able to deal with complicated geometries and boundaries. FDM is limited to rectangular shapes.\n",
      "\n",
      "* FEM quality approximation is often higher than in FDM, but this is problem dependant. (There are examples where FDM is better)\n",
      "\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Integral formulation for numerical solutions\n",
      "\n",
      "We have many integral methods that can be used to get an approximate solution to a physical problem. Let's show an example of how they work analytically.\n",
      "\n",
      "We have a simple supported beam subjected to concentrated moments at each end.\n",
      "\n",
      "<center>![image](files/resources/beam.png)</center>\n",
      "\n",
      "The governing differential equation is:\n",
      "\n",
      "$$EI \\frac{d^2y}{dx^2} - M(x) = 0$$\n",
      "\n",
      "with the boundary conditions:\n",
      "\n",
      "$$y(0)=0 \\qquad \\text{and} \\qquad y(H)=0$$\n",
      "\n",
      "Where $EI$ is the resistance of the beam to deflection and $M(x)$ is the bending moment equation, in this case has a constant value $M_o$.\n",
      "\n",
      "An approximate solution is: \n",
      "\n",
      "$$y(x)=A\\sin\\left(\\frac{\\pi x}{H}\\right)$$\n",
      "\n",
      "$A$ is a coefficient to be determined. This solution is an acceptable candidate because it satisfies the boundary conditions and also has a shape similar to the expected deflection curve. \n",
      "\n",
      "This problem has an exact solution, \n",
      "\n",
      "$$y(x)=\\frac{M_o x}{2EI}\\left(x-H\\right)$$\n",
      "\n",
      "So, let's compare the different approximations we get from the different methods.\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Variational method\n",
      "\n",
      "The variational approach involves, as is expected, the integral of a function that produces a number. The function that produces the lowest number has the property of satisfying a specific differential equation. To clarify, consider the integral:\n",
      "\n",
      "$$\\Pi=\\int_0^H \\left[\\frac{D}{2}\\left(\\frac{dy}{dx}\\right)-Qy\\right]dx$$\n",
      "\n",
      "The calculus of variations shows that the equation $y=g(x)$ which leads to the lowest numerical value for $\\Pi$, is the solution to:\n",
      "\n",
      "$$D\\frac{d^2y}{dx^2}+Q=0$$\n",
      "\n",
      "with the boundary conditions $y(0)=y_o$ and $y(H)=y_h$.\n",
      "\n",
      "Then in our example, the integral for the differential equation is:\n",
      "\n",
      "$$\\Pi=\\int_0^H \\left[\\frac{EI}{2}\\left(\\frac{dy}{dx}\\right)-M_oy\\right]dx$$\n",
      "\n",
      "Using $y(x)=A\\sin\\left(\\frac{\\pi x}{H}\\right)$ in the previous equation, integrating and making $\\partial\\Pi/\\partial A =0$ we can find the value of $A$ that minimize $\\Pi$.\n",
      "\n",
      "$$A=-\\frac{4M_oH^2}{\\pi^3EI} \\quad \\text{then,} \\quad y(x)=-\\frac{4M_oH^2}{\\pi^3EI} \\sin\\left(\\frac{\\pi x}{H}\\right)$$\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Collocation Method\n",
      "\n",
      "Impulse functions $W_i(x)=\\delta({x-X_i})$ are selected as the weighting functions (check [weighted residual method](http://nbviewer.ipython.org/github/piyueh/FEM/blob/master/notebooks/01_28_2015.ipynb)). This is equivalent to asking the integral of the residual to go to zero at specific points. The number of points selected is equal to the number of undetermined coefficients in the approximate solution.\n",
      "\n",
      "In our example the residual is obtained by substituiting the approximate solution into the PDE. \n",
      "\n",
      "$$R(x)=-EI\\frac{A\\pi^2}{H^2}\\sin\\left(\\frac{\\pi x}{H}\\right) - M_o$$\n",
      "\n",
      "due to \n",
      "\n",
      "$$\\frac{d^2y}{dx^2}=-\\frac{A\\pi^2}{H^2}\\sin\\left(\\frac{\\pi x}{H}\\right)$$\n",
      "\n",
      "There is only one undetermined coefficient. We need $R(x)=0$ at some point between $0$ and $H$. Let's choose $x=H/2$. So, doing $R(H/2)=0$ we can calculate $A$. \n",
      "\n",
      "$$A=-\\frac{M_oH^2}{\\pi^2EI} \\quad \\text{then,} \\quad y(x)=-\\frac{M_oH^2}{\\pi^2EI} \\sin\\left(\\frac{\\pi x}{H}\\right)$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Subdomain Method\n",
      "\n",
      "Each weighting function is chosen as unity, $W_i(x)=1$, over a specific region. This is equivalent to requiring the integral of the residual to vanish over an interval of the region. The number of integration intervals is equal to the number of undetermined coefficients in the approximate solution.\n",
      "\n",
      "As we mention before, in our example there is only one unknown coefficient, thus the interval must be $[0,H]$. Using the residual we calculated in the previous case, and doing $\\int R(x) dx=0$ in $[0,H]$ we get:\n",
      "\n",
      "\n",
      "$$A=-\\frac{M_oH^2}{2\\pi EI} \\quad \\text{then,} \\quad y(x)=-\\frac{M_oH^2}{2\\pi EI} \\sin\\left(\\frac{\\pi x}{H}\\right)$$\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Galerkin's Method\n",
      "\n",
      "This method uses the same functions for $W_i(x)$ that were used in the approximating equation. \n",
      "\n",
      "In our example there is only one weighting function, $W_i(x)=\\sin(\\pi x/H)$. The integral of the residual equation with the corresponding weighting function is:\n",
      "\n",
      "$$\\int_0^H \\sin\\left(\\frac{\\pi x}{H}\\right) \\left[-EI\\frac{A\\pi^2}{H^2}\\sin\\left(\\frac{\\pi x}{H}\\right) - M_o\\right]\\,dx = 0$$\n",
      "\n",
      "Integrating and solving for $A$, we get: \n",
      "\n",
      "$$A=-\\frac{4M_oH^2}{\\pi^3EI} \\quad \\text{then,} \\quad y(x)=-\\frac{4M_oH^2}{\\pi^3EI} \\sin\\left(\\frac{\\pi x}{H}\\right)$$\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Least Squares Method\n",
      "\n",
      "Utilizes the residual as the weighting function and obtains a new error term defined by:\n",
      "\n",
      "$$Er=\\int_0^H \\left[R(x)\\right]^2\\,dx$$\n",
      "\n",
      "This error is minimize the error respect to the unknown coefficients in the approximate solution. \n",
      "\n",
      "In our example, the error is:\n",
      "\n",
      "$$\\int_0^H \\left[-EI\\frac{A\\pi^2}{H^2}\\sin\\left(\\frac{\\pi x}{H}\\right) - M_o\\right]^2\\,dx$$\n",
      "\n",
      "If we integrate, minimize with respect to $A$ ($\\partial Er/ \\partial A = 0$) and then solve for $A$, we get:\n",
      "\n",
      "\n",
      "$$A=-\\frac{4M_oH^2}{\\pi^3EI} \\quad \\text{then,} \\quad y(x)=-\\frac{4M_oH^2}{\\pi^3EI} \\sin\\left(\\frac{\\pi x}{H}\\right)$$\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### About the percent error\n",
      "\n",
      "It is not possible to say which method is the more accurate because the error depends on the approximating function and the equation being solved. For our example, the percent error in deflection is:\n",
      "\n",
      "<center>![image](files/resources/error_beam.png)</center>\n",
      "\n",
      "It seems that Variational, Galerkin's and Least squares methods are more accurate than the others. However, it is possible to find a collocation point which produces a maximum deflection that agrees with the exact value. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "*The important message is that the numerical solution of a differential equation can be formulated in terms of and integral. The integral formulation is a basic characteristic of the finite element method.*"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###References:\n",
      "\n",
      "* Applied finite element analysis. *Larry J. Segerlind*. Second edition (1984)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For code examples (Python), and a really good tutorial on FEM go to:\n",
      "\n",
      "[Introduction to finite element methods](http://hplgit.github.io/INF5620/doc/pub/main_fem.pdf)\n",
      "It has pieces of code in the pdf, but the codes are also available on github.\n",
      "[codes](https://github.com/hplgit/num-methods-for-PDEs)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": []
    }
   ],
   "metadata": {}
  }
 ]
}